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In most theoretical descriptions of collective strong coupling of organic molecules to a cavity 
mode, the molecules are modeled as simple two-level systems. This picture fails to describe the rich 
structure provided by their internal rovibrational (nuclear) degrees of freedom. We investigate a 
first-principles model that fully takes into account both electronic and nuclear degrees of freedom, 
allowing an exploration of the phenomenon of strong coupling from an entirely new perspective. First, 
we demonstrate the limitations of applicability of the Born-Oppenheimer approximation in strongly 
coupled molecule-cavity structures. For the case of two molecules, we also show how dark states, 
which within the two-level picture are effectively decoupled from the cavity, are indeed affected by 
the formation of collective strong coupling. Finally, we discuss ground-state modifications in the 
ultra-strong coupling regime and show that some molecular observables are affected by the collective 
coupling strength, while others only depend on the single-molecule coupling constant. 

PACS numbers: 71.36.-fc, 78.66.Qn, 82.20.Kh, 73.20.Mf, 


I. INTRODUCTION 

Strong coupling in quantum electrodynamics is a well- 
known phenomenon that occurs when the coherent energy 
exchange between a light mode and quantum emitters 
is faster than the decay and decoherence of either con¬ 
stituent [1, 2]. The excitations of the system are then 
hybrid light-matter excitations, so-called polaritons, that 
combine the properties of both constituents. Exploiting 
these properties enables new applications such as polariton 
condensation under collective strong coupling to excitons 
(excited electron-hole pairs) in semiconductors [3, 4] and 
organic materials [5-7]. Organic materials present a par¬ 
ticularly favorable case, as the Frenkel excitons in these 
materials possess large binding energies, large dipole mo¬ 
ments, and can reach high densities. This enables Rabi 
splittings Ofl (the energy splitting between the polaritons) 
up to more than 1 eV [8-10] , a significant fraction of the 
uncoupled transition energy. These properties allow for 
strong coupling to many kinds of electromagnetic (EM) 
modes [11], such as cavity photons [8, 9, 12], surface plas- 
mon polaritons [13-16], surface lattice resonances [17, 18], 
or localized surface plasmons [19, 20] . 

While organic molecules are thus uniquely suited to 
achieving strong coupling, they are not simple two-level 
quantum emitters, but rather have a complicated level 
structure including not only electronic excitations, but 
also rovibrational degrees of freedom (schematically de¬ 
picted in Eig. 1). It has been experimentally demon¬ 
strated that strong coupling can modify this structure, in 
the sense that material properties and chemical reaction 
rates change [21-23]. However, the models used to de¬ 
scribe strong coupling are often focused on macroscopic 
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FIG. 1. Illustration of energy level structure of a bare com¬ 
plex molecule and the hybrid states that result in the strong 
coupling regime with a photonic mode of energy hoJc, resonant 
with the molecular excitation. 


descriptions [24], and most microscopic models do treat 
organic molecules as two-level systems (see [25] for a re¬ 
cent review). When the rovibrational degrees of freedom 
are taken into account, this is often done using effective 
decay and dephasing rates [26] , with a few works explic¬ 
itly including a phononic degree of freedom [27, 28]. All 
of these approaches only provide limited insight into the 
effects of strong coupling on molecular structure. 

In the present work, we thus aim for a microscopic 
description of strong coupling with organic molecules. 
We introduce a simple first-principles model that fully 
describes nuclear, electronic and photonic degrees of free¬ 
dom, but can be solved without approximations. This 
allows us to provide a simple picture for understanding 
the induced modification of molecular structure. 

In section H, after introducing the model, we discuss 
under which conditions and in which form the Born- 
Oppenheimer approximation (BOA) [29, 30] is valid in 
the strong coupling regime for a single molecule. The 
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BOA is widely used in molecular and solid state physics 
and quantum chemistry, and provides a simple picture of 
nuclei moving on effective potential energy surfaces (PES) 
generated by the electrons, which underlies most of the 
current understanding of chemical reactions [30]. How¬ 
ever, the BOA depends on the separation of electronic and 
nuclear energy scales, i.e., the fact that electrons typically 
move much faster than nuclei. It could thus conceivably 
break down when an additional, intermediate timescale is 
introduced under strong coupling to an EM mode. The 
speed of energy exchange between field and molecules is 
determined by the Rabi frequency Hr, and typical exper¬ 
imental values of hundreds of meV land squarely between 
typical nuclear (~ 100 meV) and electronic (~ 2 eV) en¬ 
ergies. We show that the BOA indeed breaks down at 
intermediate Rabi splittings, but remains valid when CIr 
becomes large enough. For cases where it breaks down, 
we show that the non-BO coupling terms can be obtained 
to a good approximation without requiring knowledge of 
the electronic wave functions. 

In section III, we focus on the effects of strong cou¬ 
pling when more than one molecule is involved, using 
two molecules as the simplest test case. In experiments, 
strong coupling is achieved by collective coupling to a large 
number of molecules, under which the Rabi frequency is 
enhanced by a factor of \/N. The number of emitters N 
is on the order of > 10^ within cavities [8-10, 12], with 
plasmonic nanoparticle modes allowing to reduce this to 
N ~ 100 [20]. In this context, it is well known that only 
a small fraction of the collective electronic excitations are 
strongly coupled ]25, 31, 32], with a large number of “dark” 
or “uncoupled” modes that show no mixing with the EM 
mode and no energy shift. We show that even these dark 
modes are affected by strong coupling, with the nuclear 
motion of separated molecules becoming correlated. 

In section IV, we focus on the so-called ultrastrong 
coupling regime, where the Rabi frequency reaches a sig¬ 
nificant fraction of the electronic transition energy, as 
achieved in experiments. In this regime, not just excited- 
state, but also ground-state properties are modified—for 
example, the ground state acquires a photonic contribu¬ 
tion [24, 33]. Accordingly, we discuss whether ground- 
state chemical properties of organic molecules could be 
modified by strong coupling. This also allows us to par¬ 
tially answer the open question what strong coupling 
means for modifications of chemical structure [34], i.e., 
whether “all” molecules are modified by it, or only a small 
subset, or whether we necessarily have to invoke collective 
modes even when discussing “single-molecule” effects. We 
show that while some observables, such as energy shifts, 
are determined by the collective Rabi frequency, but other 
observables, such as the shift in ground-state bond length, 
are instead determined by the single-molecule coupling 
strength oc Qr/'/N- 

For simplicity, we only treat a single EM mode and 
completely neglect dissipation in the following. We use 
atomic units unless stated otherwise (47reo = h = = 

e = 1, with electron mass me and elementary charge e). 


II. SINGLE MOLECULE 

In this section, we introduce our model for a single 
molecule coupled to an EM mode. Due to the exponen¬ 
tial scaling with the degrees of freedom, solving the full 
time-independent Schrodinger equation for an organic 
molecule without the BOA is an extremely challenging 
task that even modern supercomputers can only handle 
for very small molecules. We thus employ a reduced- 
dimensionality model that we can easily solve, both for 
the bare molecule and after coupling to an EM mode. 


A. Method 

1. Bare molecule 


We work within the single-active-electron approxima¬ 
tion (SAE), in which all but one electron are frozen around 
the nuclei, and additionally restrict the motion of the ac¬ 
tive electron to one dimension, x. Furthermore, we only 
treat one nuclear degree of freedom, the reaction coor¬ 
dinate R. This could correspond to the movement of a 
single bond in a molecule, but can equally well represent 
collective motion, e.g., the breathing mode of a carbon 
ring. The effective molecular Hamiltonian then highly 
resembles that of a one-dimensional diatomic molecule, 

Hm = + Te + yen{x]R) -f (1) 

where T„ = and Tj. = 4“ nuclear and electronic 

kinetic energy operators (with P, p the corresponding 
momenta), and M is the nuclear mass. The potentials 
Venix,R) and Vnn{R) represent the effective electron- 
nuclei and internuclear interactions, where we assume two 
nuclei located at x = ±Rl2. These potentials encode 
the information about the frozen electrons as well as the 
nuclear structure of the molecule, and can be adjusted to 
approximately represent different molecules. 

The electron-nucleus interaction Ven contains the inter¬ 
action of the active electron with each nucleus, as well 
as with the frozen electrons surrounding it. Assuming 
a nuclear charge of Z, we have 2Z — 1 frozen electrons 
distributed across the two nuclei. For large distances, the 
active electron should thus feel a Coulomb potential with 
an effective charge of | from each nucleus. Conversely, at 
very small distances, the active electron is not affected by 
the cloud of frozen electrons and feels an effective charge 
of Z. Since we are working within one dimension, we use 
a soft Coulomb potential to take into account that the 
electron avoids the singularity at the nucleus. We choose 
a simple model potential fulfilling these conditions: 


Ven{r) 


1 

2 


+ {Z-l)e-^ 


( 2 ) 


where a is the softening parameter, rg describes the lo¬ 
calization of the frozen electrons around the nucleus, and 
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r is the electron-nucleus distance. The total potential is 
thus Ven{x, R) = Ve„{\x - R/2\) + Ven{\x + R/2\). 

The internuclear potential VnniR) represents the inter¬ 
action between the nuclei and the 2Z — 1 frozen electrons, 
i.e., the ground state potential energy surface of the molec¬ 
ular ion. We model this surface by a Morse potential 

VnniR) = De (l - , (3) 

which adds three new parameters: the dissociation energy 
De, the equilibrium distance Rq and the width of the 
potential well A. By tuning the seven free parameters we 
have at our disposal (M, Z, a, tq, Df., Rq and A), we 
can approximately fit both the electronic and vibrational 
structure and absorption spectrum to those of real organic 
molecules. 

We can now solve the stationary Schrodinger equation 
Hm'^iXjR) = E'i>ix,R) for the bare-molecule Hamilto¬ 
nian Eq. (1) without further approximations by repre¬ 
senting Hm on a two-dimensional grid in x and R. For a 
bare molecule, the results are virtually identical to those 
obtained within the BOA and thus not shown here. 

We next give a short description of the Born- 
Oppenheimer approximation for completeness (see [29, 30] 
for more details). As mentioned above, the basic idea is 
to exploit the separation between nuclear and electronic 
timescales and to assume that the electrons perfectly 
follow nuclear rearrangements without changing state 
(i.e., adiabatically). This is achieved by separating the 
Hamiltonian into the nuclear kinetic energy T„ and an 
electronic Hamiltonian Heix] R) = Hmix,R) — Tn that 
only depends on R parametrically. Diagonalizing He 
yields a set {(j>k} of electronic eigenstates for every R, 
with Heix] R)4>kix', R) = Ek{R)4>kix\ R)- Without loss of 
generality, each total eigenstate T* can be represented by 
4'*(a;,i?) = '^k^kix\R)XkiR)- Inserting this expansion 
into the Hamiltonian Eq. (1) leads to a set of coupled 
differential equations, 

ifn + Ek)xUR) + E = ^xUR), (4) 

k' 

with nuclear motion taking place on potential energy 
surfaces (PES) Ek{R) that are coupled through correc¬ 
tion terms = {(t>k\fn\(t)k')x + {(t>k\^\(l}k')xP, where 
the subscript x indicates that the integration in the 
brakets is only over the electronic coordinate. The 
Born-Oppenheimer approximation now consists in ne¬ 
glecting the intersurface couplings (7^^ , which can be 
shown to be small when the electronic levels are well- 
separated. This gives a set of independent PES Ek{K) 
on which the nuclei move, where each eigenstate is a 
product of a single electronic and nuclear wave function, 
'^\ix,R) = 4>kix] R)XkiR)■ Ths different nuclear func¬ 
tions on each electronic curve correspond to rotational 
or vibrational excitation. This picture of nuclear motion 
on PES is extremely powerful and underlies most of the 
current understanding of chemical reactions [30]. The 
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FIG. 2. Bare-molecule potential energy surfaces of the two 
first electronic states in the BOA for (a) the rhodamine 6G-like 
model molecule and (c) the anthracene-like model molecule. 
The vibrational levels and associated nuclear probability densi¬ 
ties are represented on top of the PES. (b) and (d): Absorption 
spectrum for the (b) R6G-like and (d) anthracene-like molecule 
in arbitrary units. 


question of its validity in the strong coupling regime is 
thus of central importance for the possible modification of 
chemical reactions and structure through strong coupling. 

In the following, we will focus on two model molecules, 
which approximately reproduce the absorption spectra 
of rhodamine 6G (R6G) and anthracene molecules that 
are commonly used in experimental realizations of strong 
coupling [12, 15, 17]. Only the first two PES, correspond¬ 
ing to the ground EgiR) and first electronically excited 
Ee{R) state, play a role in the results discussed in the 
following. They are shown in Fig. 2a and Fig. 2c, to¬ 
gether with the nuclear probability densities |x(i?)p for 
the lowest vibrational levels on each PES. Importantly, 
the two models differ significantly in two relevant quanti¬ 
ties: the vibrational mode frequency ujyib and the offset 
AR, i.e., the change in equilibrium distance between the 
ground and excited PES. This offset is related to the 
strength of the electron-phonon interaction and influences 
the Stokes shift between emission and absorption [35] . 
The R6G-like model has relatively small vibrational spac¬ 
ing uJvib « 70 meV and small offset AR « 0.018 a.u., while 
the anthracene-like model has large vibrational spacing 
UJvib ~ 180 meV and large offset AR « 0.092 a.u.. Ac¬ 
cordingly, their absorption spectra (Fig. 2b and Fig. 2d, 
obtained from Eq. (6) below) are qualitatively different, 
with anthracene showing a broader absorption peak with 
well-resolved vibronic subpeaks. 
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FIG. 3. Strongly coupled electronic PES (solid lines) in the 
singly excited subspace, for the anthracene-like molecule for 
(a) g = 0.001 a.u. and (b) g = 0.008 a.u.. The dashed lines 
show the corresponding uncoupled states: A molecule in the 
first excited state, Ee{R), and a molecule in the ground state 
with one photon present, Eg{R) + ujc- 


2. Molecule-photon coupling 


We now add a photonic mode and its coupling to the 
molecule (within the dipole approximation) into the molec¬ 
ular Hamiltonian, 

Hmc = Hm+uica^a ++a), (5) 


where /t is the dipole operator of the molecule {fi = x va. 
our case), and a are the creation and annihilation op¬ 
erators for the bosonic EM held mode, Wc is its frequency, 
and g is the coupling strength constant, given by the elec¬ 
tric held amplitude (along the x-axis) of a single photon. 
In the following, we always set the photon energy Wc to 
achieve “zero detuning”, with Wc at the absorption maxi¬ 
mum of the molecule. This gives ojc « Ee{Re) — Eg{Re), 
where Re is the equilibrium position at which Eg[R) has 
its minimum. 

To provide some context for the held strengths used 
in the following, we note that for a typical microcavity 


hujc 

2 ^ 


with a mode volume E « A^, one obtains 5 = y 

1.34 X a.u. (for ujc given in eV). However, for an 

effective mode volume close to the current record achieved 
in plasmonic nano-antennas, V « 1.3 x 10“^Ag [36], the 
single-particle coupling reaches g « 3.72 x 10“^a;^ a.u. 
(wc again in eV). We furthermore note that the ground- 
to-excited state dipole transition moments of our model 
molecules are on the order of 1 a.u. « 2.54 D, i.e., almost 
an order of magnitude smaller than in typical organic 
molecules [37]. 

Compared to the bare-molecule case, the Hamiltonian 
now includes a new degree of freedom, the photon number 
n € { 0 , 1 , 2 ,...}, with the system eigenstates dehned 
by R) = E^{x,n, R). As discussed above, 

the typical energies associated with strong coupling in 
organic molecules are somewhere between the nuclear and 
electronic energies. A priori, this suggests two options of 
performing the BOA: the additional terms introduced by 
the photonic degree of freedom could be grouped either 
with the “slow” nuclear motion or with the “fast” electronic 


Hamiltonian. However, as the photon couples to the 
electron, grouping it with the nuclear terms necessarily 
leads to additional off-diagonal terms in Eq. (4), and 
no independent PES on which the nuclei move could be 
obtained. Consequently, the only way to maintain the 
usefulness of the BOA is to include the photonic degree 
of freedom within the electronic Hamiltonian, leading to 
a new set of “strongly coupled PES”. 

We first focus on the singly excited subspace, within 
which the splitting between polaritons is observed. Here, 
either the molecule is electronically excited and no pho¬ 
tons are present, or the molecule is in its electronic ground 
state and the photon mode is singly occupied. At zero 
coupling (g — 0), this gives two uncoupled PES {Ee{R) 
and Eg{R) -l-Wc, dashed curves in Fig. 3) that cross close 
to Re for our choice of Wc. When the electron-photon 
coupling is non-zero but small, a narrow avoided crossing 
develops instead (solid lines in Fig. 3a), while for large 
coupling strengths, the energy exchange between pho¬ 
tonic and electronic degrees of freedom is so fast that we 
observe two entirely new PES (Fig. 3b), which can not 
be easily associated with either of the uncoupled PES. 
Instead, they become hybrid polaritonic PES that con¬ 
tain a mixture of electronic and photonic excitation, the 
hallmark of the strong coupling regime. 

As discussed above, the BOA is known to be valid when 
two PES are sufficiently separated from each other. This 
implies that the BOA breaks down when g is small and 
the two PES possess a narrow avoided crossing. This 
in itself is not a surprising result—when the electron- 
photon coupling is very small, the system is not even 
in the strong coupling regime, and the photon mode 
is better treated as a small perturbation. Fortunately, 
the weak coupling regime is also not interesting from 
the standpoint of understanding or modifying molecular 
structure through strong coupling. The real question thus 
must be: How strong does the electron-photon coupling 
have to be for the BOA to be valid, and is this condition 
fulfilled for realistic experimental parameters? In order to 
better quantify the agreement between the BOA and the 
full solution, we next turn to an easily measured physical 
observable: the absorption spectrum. 


B. Absorption 


The absorption cross section at frequency uj can be 
calculated from the polarizability as [38, 39] 


<7{UJ) 



l(V’fclAIV^o)P 

ujk — ojQ — u: — 


( 6 ) 


where the sum runs over all eigenstates \'4>k) of the system 
(with energies uik) and |'0o) is the ground state. As we do 
not include incoherent processes in our calculation, this 
would give d-like peaks in the absorption cross section. 
In the plots shown in the following, we instead introduce 
a phenomonological width representing losses and pure 
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FIG. 4. Absorption cross sections of a single molecule, calcu¬ 
lated using the full Hamiltonian without approximation (solid 
green lines) and within the BOA (dashed black lines). Results 
are shown for the (a) R6G-like and (b) anthracene-like model 
molecules, for several values of the coupling strength g. 


dephasing by setting e to a small non-zero value, such that 
the absorption cross section becomes a sum of Lorentzians. 
For the bare-molecule case without coupling to an EM 
mode, the absorption spectra of our two model molecules 
approximately agree with those of R6G (Fig. 2b, [17]) 
and anthracene (Fig. 2d, [12]). 

In Fig. 4, we compare the absorption cross sections 
under strong coupling as obtained from a full calculation 
without approximations to those obtained within the 
BOA, for a range of coupling strengths g to the EM 
mode. Even for relatively small g, the BOA is found to 
agree almost perfectly with the full results for the R6G- 
like molecule with small vibrational spacing (Fig. 4a). 
However, for the anthracene-like molecule with a high- 
frequency vibrational mode and large offset zAR, the BOA 
only agrees with the full resul for relatively large values 
of g, where the Rabi splitting flu (as defined by the 
energy difference between the two “polariton” peaks in 
the absorption spectrum) is appreciably larger than the 
vibrational frequency « 180 meV (Fig. 4b). 

This qualitative observation can be quantified by com¬ 
paring the correction terms (7^^ in Eq. (4) with the 
energy difference between the anticrossing PES at the 
point of closest approach. In appendix A, we present 
a model that achieves this without any explicit knowl¬ 
edge of the electronic wave functions. It relies on the 
observation that close to the anticrossing, the coupled 
(polariton) states switch character between the two un¬ 
coupled states, while the “intrinsic” R-dependence of the 
uncoupled electronic states can be neglected. The cor¬ 
rection terms (7^^ can then be obtained just from the 
knowledge of Eg{R), Ee(R), and g,gg{R), where geg{R) is 
the electronic transition dipole between the ground and 
excited state. By approximating Eg{R) and Ee{R) as 
harmonic oscillators, the correction terms can be analyti¬ 
cally evaluated and are found to be negligible under the 
condition that ARw^j{,/0|j is small compared to the nu¬ 


clear momentum of the relevant eigenstates. This demon¬ 
strates that the model molecules present two opposite 
cases for the applicability of the BO approximation: our 
R6G-like molecule has a relatively small vibrational spac¬ 
ing uj-uib ~ 70 meV and small electron-phonon coupling, 
AR « 0.018 a.u., while our anthracene-like model molecule 
has a large vibrational spacing ojyib « 180 meV and large 
electron-phonon coupling, AR « 0.092 a.u.. We note 
that in many experiments involving organic molecules, 
ilfi > 500 meV [9, 10] is significantly larger than typical 
vibrational frequencies ujyib ^ 200 meV [40] . This shows 
that the intuitive picture of nuclear dynamics unfolding 
on uncoupled Born-Oppenheimer potential energy sur¬ 
faces can often be applied to understand the modification 
of molecular chemistry induced by strong coupling. Ad¬ 
ditionally, even when the BOA breaks down, the model 
presented in appendix A can be used to obtain the non- 
BO coupling terms without requiring knowledge of the 
electronic wave functions. The only necessary input are 
the uncoupled PES and the associated transition dipole 
moments. Even for relatively large molecules, these can 
be obtained using the standard methods of quantum 
chemistry or density functional theory. 


III. TWO MOLECULES 

In the previous section, we showed that on the single¬ 
molecule level, the BOA is valid as long as the Rabi 
splitting rin is large enough. However, current experi¬ 
ments are performed with large numbers of molecules, 
where coherent superpositions of electronic excitations 
(bright “Dicke states” [41]) couple strongly to the pho¬ 
tonic mode(s), while other superpositions give uncoupled 
or “dark” modes. It is thus important to consider if and 
how our conclusions have to be modified when more than 
a single molecule is involved in strong coupling. 

For later reference, we give a quick overview of the 
theory when using an ensemble of two-level emitters 
coupled to a photonic mode, i.e., the many-particle 
Jaynes-Gummings (JG) model [42], also known as the 
the Tavis-Gummings model [43]. Its Hamiltonian within 
the rotating-wave approximation is 

Hjc = Wed’^'a + ^ Wic|ci + '^gi{acl + a^c*), (7) 

i i 

where u>i is the energy of emitter i with destruction 
(creation) operator Ci (c|^), and the gi describe the 
emitter-photon couplings. For identical emitters (wi = 
= g)y the resulting eigenstates in the single¬ 
excitation subspace are given by two polaritons |±) = 
;^(dt|0) ± \B)), symmetric and antisymmetric combina¬ 
tions of the photonic mode with the emitter bright state 
1^) = ^LjC-lO)- At zero detuning (w^ = w^), the 
polariton energies are given by u!± = ujm ± ^lnl2, where 
Vlji = 2g\/N is the collective Rabi splitting. The N — 1 
superpositions of emitter states orthogonal to |R) are dark 
















6 



FIG. 5. (a) Uncoupled potential energy surfaces of two anthracene-like molecules in the singly excited subspace: Esgo{Ri, R 2 ) 
(orange), Eego{Ri, R 2 ) (blue), and Eggi{Ri, R 2 ) (green), (b) Coupled PES for g = 0.002 a.u. and (c) g = 0.013 a.u., corresponding 
to the lower polariton (orange), dark state (blue), and upper polariton (green). For clarity, only parts where Ri < R 2 are shown 
(note that the system is symmetric under the exchange i?i -f-)- i? 2 ). 


states that are not coupled to the photonic mode, with 
energies identical to the uncoupled emitters, ujds = 
Note that in configurations with many photonic modes 
(e.g., planar cavities), more than one emitter state is 
coupled to the photonic mode (typically at low in-plane 
momentum), but there remain many uncoupled (dark) 
modes at higher in-plane momentum [25, 32]. There is 
an ongoing discussion in the literature on whether the 
“dark” modes are affected by strong coupling as well, or 
whether they should be thought of as completely unmodi¬ 
fied emitter states. We will show below that when taking 
the internal structure of the emitters (molecules) into 
account, even the “dark” modes are affected by strong 
coupling and the nuclear dynamics of separate molecules 
become correlated. 

We now treat the case of two model molecules, which 
can still be solved exactly within our approach, but which 
displays many of the effects of many-molecule strong 
coupling. As in the JC model, we assume that the two 
molecules both couple to the same photonic mode, but 
do not directly interact with each other, giving 

= WcO^'a -b { h ^'> + + a)) , (8) 

i=i.2 

where the superscripts j indicate the molecule on which 
the operator acts. Directly diagonalizing this Hamilto¬ 
nian in the “raw” basis {xi, i?i, 0 : 2 , i? 2 , is already a 
formidable computational task for typical grid sizes. We 
thus calculate the full solution by first diagonalizing the 
single-molecule Hamiltonian, 

including only a relevant subset of eigenstates {k} for 
each molecule in the total basis {ki,k 2 ,n}. The number 
of necessary eigenstates to obtain completely converged 
results is quite small (« 30 per molecule). However, this 
approach only provides limited insight into the dynam¬ 


ics of the strongly coupled system, especially regarding 
nuclear motion. 

We thus again apply the Born-Oppenheimer approxi¬ 
mation by separating the nuclear kinetic energy terms and 
diagonalizing the remaining Hamiltonian parametrically 
as a function of i?i and i? 2 - Similar to above, instead 
of working in the {xi,X 2 ,n\ basis for each combination 
(i?i,i? 2 ), we prediagonalize the single-molecule electronic 
Hamiltonian Hf.{x\R) = '^^Ek{R)\k{R)){k{R)\, where 
(for the cases discussed here) the sum only has to include 
the ground and first excited states to achieve convergence, 
k S {g, e}. If we additionally allow at most one photon in 
the system, n G {0,1}, we obtain an 8 x 8 Hamiltonian 
for each combination of nuclear coordinates i?i, i? 2 . 

The electronic Hamiltonian consists of all possible com¬ 
binations of electronic states Eg, E^ of the two molecules 
with 0 or 1 photons. A further simplification is achieved 
by taking into account that the Hamiltonian conserves 
parity H = (—l)’^i+'^ 2 +"j with TTj the parity of the state 
of molecule j (gerade or ungerade). For large coupling g, 
this separation by parity avoids some accidental degenera¬ 
cies between uncoupled PES and thus improves the BOA. 
We now obtain two independent 4x4 Hamiltonians, 


.Heven (.^1, R 2 ) 


^odd(-Rl, .R2) 


/^ggo 


gd(^) 



Eegl 

0 


gS‘^'> 

0 

Egel 

gd^^^ 

\ 0 



EeeO / 

/ddggl 

gS^'> 



gS^'> 

^egO 

0 



0 

EgeO 


\ 0 

gS^'> 


Eeel ) 


(9a) 


(9b) 
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R6G-Ukc molecule 



FIG. 6. Absorption cross section of two molecules driven 
coherently, calculated using the full Hamiltonian without ap¬ 
proximation (solid green lines) and within the BOA (dashed 
black lines). Results are shown for the (a) R6G-like and (b) 
anthracene-like model molecules, for several values of the cou¬ 
pling strength g. The values of g are scaled by l/\/2 with 
respect to the single-molecule case (Fig. 4) in order to obtain 
the same total Rabi frequency flij. 


where the uncoupled PES are represented by the com¬ 
pact notation Eabn = Ea{Ri) + Eb{R 2 ) + nWc, and the 
single-molecule dipole transition moment between the 
ground and first excited state is denoted by = 
{(j)g{Rj)\jl\(j)e{Rj)) ■ Diagonalizing these Hamiltonians for 
each (i?i,i? 2 ) results in a set of strongly coupled two- 
dimensional PES. In Eig. 5, we show the three surfaces 
in the single-excitation subspace, corresponding to the 
three lowest states of Eq. (9b). For zero molecule-photon 
coupling {g = 0, Fig. 5a), there are now a number of 
one-dimensional seams where the three PES cross. When 
the molecule-photon coupling is turned on, these crossings 
again turn into avoided crossings, as shown in panels (b) 
and (c) for two different coupling strengths g. Follow¬ 
ing the conventions used in the Jaynes-Cummings model, 
we label the three coupled PES in order of energy as 
the “lower polariton PES”, the “dark-state PES”, and the 
“upper polariton PES”. 

We first address the applicability of the BOA, which 
breaks down when two PES are close in energy, for the case 
of two molecules. Within the single-excitation subspace 
(which determines the linear properties of the system, 
such as absorption), there are now a range of (avoided) 
crossings. They occur when i) all three surfaces approach 
each other Eggi « Egeo ~ Ef,gQ, ii) the photonically ex¬ 
cited PES is close to only one of the electronically excited 
PES, Eggi « Egeo Or Eggi « Eggo, Or hi) only the two 
electronically excited states cross, Eg^o « Eego- Case i) 
corresponds to the JC model at zero detuning, giving 
the two polaritonic PES at energy shifts of and 

an additional dark state that is unshifted from the bare- 
molecule case. The BOA in this region is thus valid for 
similar conditions as in the single-molecule case, although 
the PES separation is reduced by half due to the addi¬ 


tional dark-state surface. Case ii) corresponds exactly to 
the single-molecule case, with the second molecule act¬ 
ing as a “spectator” that only induces additional energy 
shifts. The BOA should thus again be valid for similar 
conditions as with a single-molecule, albeit with the cou¬ 
pling reduced by l/-\/2 for a fixed total Rabi splitting. 
Finally, case iii) presents the biggest challenge, as the 
two electronically excited PES, E^go and Eg^o, are not di¬ 
rectly coupled, but only split indirectly through coupling 
to the photonically excited surface Eggi. The splitting 
between the two surfaces is thus small for large detuning, 
AE « (gd)^ /4:{Eggi — Eego). This is clearly observed in 
Fig. 5b along the line i?i = R 2 , where the dark state PES 
almost touches the upper polariton PES for small Rs and 
the lower polariton PES for large Rs. 

The discussion above implies that for almost any cou¬ 
pling strength, there will be regions in the nuclear con¬ 
figuration space where the BOA breaks down. 

However, not all parts of the PES are visited by the nu¬ 
clei during a given physical process. To explicitly check 
the BOA in the subspace relevant for polaritonic physics, 
in Fig. 6 we thus again compare the absorption with that 
obtained by a full diagonalization of the Hamiltonian 
Eq. (8). Compared to the single-molecule case, many 
more molecular levels are present in the system, leading 
to small changes in the absorption spectra compared to 
the single-molecule case. In order to properly compare 
the results, we take into account the '/N scaling of the 
total Rabi frequency and reduce the coupling strengths 
by to produce the same total splitting. The BOA is 
shown to again become valid for large enough coupling, 
but the minimum coupling required is increased compared 
to that for a single molecule. In the common case of slow 
nuclear motion, as for our R6G-like model in Fig. 6a, the 
BOA already is valid for relatively small Rabi splitting 
of « 250 eV. However, in the anthracene-like case 
of very fast vibrational motion. Fig. 6b, the BOA still 
does not give perfect agreement with the full model for 
g = 0.0057 a.u. (O^ « 600 meV), and agreement is only 
reached at roughly twice that value. 

Having established the validity of the BOA for many 
relevant cases and Rabi splittings comparable to experi¬ 
mental values, we now discuss the implications of collec¬ 
tive strong coupling for the internal molecular (nuclear) 
dynamics. Note that this question can by design not 
be addressed within the JC model, where emitters are 
two-level systems without any internal degrees of free¬ 
dom. In contrast, the BOA provides a straightforward 
approach to this problem. Any two-dimensional PES can 
be decomposed into a sum of independent single-molecule 
potentials, plus a remainder that describes the coupling 
between the nuclear motion of the molecules, 

E(i?i,i?2) =Ei(i?i)+E2(R2)+Ei2(i?l,i?2). (10) 

The nuclear motion of two molecules is independent 
if and only if the coupled part Ei 2 {Ri, R 2 ) is identically 
zero. In order to quantify this coupling, we expand each of 
the coupled PES in the single-excitation subspace around 


















FIG. 7. Coupling between nuclear motion in different 
molecules for the lower (LP) and upper polariton (UP) and 
dark-state (DS) PES. Results are shown as the ratio /3 /q; be¬ 
tween the prefactors of the offdiagonal R 1 R 2 and diagonal R^ 
terms in Eq. (11), for the R6G-like model molecule. 

its minimum giving 

E{Ri,R2) « Eo + a5Rl + a5Rl -k (H) 

with Eq = E{R\,R 2 ) and 5Ri = Ri — R^. Note that due 
to symmetry under the exchange iii o i? 2 , the prefactor 
a is the same for 5R\ and SR^- As can be seen in Fig. 7, 
both the polariton and even the dark state PES show sig¬ 
nificant coupling of the nuclear degrees of freedom, with 
values of /3/a on the order of a few percent for values of 
9 ^ 0.01 a.u. giving realistic Rabi splittings of ;< 1 eV 
(see Fig. 6). Interestingly, the coupling is much larger 
for the lower polariton state than for either the upper 
polariton or the dark state, and decreases with increasing 
g for all three PES. We therefore conclude that even dark 
states that have negligible mixing with photonic modes 
are affected by strong coupling, in the sense that the 
nuclear degrees of freedom of separate molecules behave 
like coupled harmonic oscillators, and their motion be¬ 
comes correlated. This implies that, e.g., local excitation 
of nuclear motion within one molecule could affect the 
nuclear motion in another, spatially separated molecule, 
even when no photon is ever present in the system. Note 
that these results apply within the singly-excited sub¬ 
space, i.e., the coupled nuclear motion is only observed 
when electronic excitation is present, not in the ground 
state. In the next section, we discuss which modifications 
of the ground state properties could be observed in the 
ultrastrong coupling regime. 

IV. ULTRASTRONG COUPLING AND 
GROUND STATE MODIFICATIONS 

Up to now, we focused on the molecular properties in 
the singly excited subspace, which are probed in linear re¬ 
sponse measurements such as absorption and transmission. 


and where the effect of strong coupling is immediately ap¬ 
parent. However, when the Rabi frequency, i.e., the energy 
exchange rate between the molecules and the photonic 
mode, becomes significant compared to the frequencies 
of these two modes, the so-called ultra-strong coupling 
regime is entered [24, 33]. In this regime, the rotating 
wave approximation for the emitter-light interaction (un¬ 
der which the number of excitations is conserved) becomes 
invalid. In our approach, the rotating wave approxima¬ 
tion is not used, and we can thus naturally explore the 
ultrastrong coupling regime. One of its most intriguing 
predictions is that even the ground-state properties of the 
system should be significantly modified. Eor example, the 
ground state is shifted in energy and acquires a photonic 
component, such that a number of virtual photons are 
present in the system even without any external excita¬ 
tions. This raises the question of how the internal degrees 
of freedom of organic molecules are affected when this 
regime is entered. 

The BOA is well-suited to explore this regime. In con¬ 
trast to the singly-excited subspace, where narrow avoided 
crossings can affect its validity, the ground-state PES is 
energetically well-separated from all other PES. This 
remains true even under ultrastrong coupling, and conse¬ 
quently the BOA is valid for all coupling strengths. The 
ground state potential energy surface Eg (R) is coupled to 
the doubly excited surface Ee{R) + 00 ^ (cf. Eq. (9a)), with 
the strongly coupled ground state PES given to lowest 

order by E^^iR) « EgiR) - + ^(9^^ 

Ground state properties such as the bond length are 
determined by the shape of the PES. The largest mod¬ 
ification can thus be expected when the i?-dependence 
of the ground and excited PES is as different as possible. 
This occurs for large electron-phonon coupling, i.e., a 
large value of AR, such as in our anthracene-like molecule. 
For a coupling strength oi g = 0.016 a.u., correspond¬ 
ing to a Rabi splitting of « 1.2 eV in absorption, 
we obtain a shift in the ground state bond length of 
Ai?o ~ 0.84 mA= 84 fm. While small, such a change 
in bond length could be detectable using X-ray absorp¬ 
tion fine structure spectroscopy or X-ray crystallography, 
which can obtain few- or even sub-femtometer precision 
for measuring bond length shifts [44, 45] . 

However, the previous paragraph only applies for a 
single molecule under strong coupling. Repeating the cal¬ 
culation using two molecules and taking into account the 
reduced single-molecule coupling strength (for fixed Rabi 
splitting, Qn oc VNg) reveals a reduction of the bond- 
length change by a factor of two, ARq™"* « 0.42 mA. 
This is confirmed by using a similar analytical model 
as presented in appendix A, in which the bare-molecule 
potential energy surfaces are approximated as harmonic 
oscillators. Due to the simple analytical structure, pertur¬ 
bation theory can be applied to obtain results for arbitrary 
numbers of molecules, and the change in ground-state 
bond length is found to be proportional to the squared 
single-molecule coupling g^, not to the squared collective 
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Rabi frequency r2|.. We note that in contrast, the ground- 
state energy shift is indeed determined by the collective 
coupling strength, AEq oc In realistic experimental 
configurations involving large numbers of molecules, the 
change in ground state bond length is thus expected to 
be minuscule and extremely challenging to measure. This 
demonstrates that the influence of strong coupling on 
any specific observable is not immediately obvious, and 
has to be checked case by case. For some properties, the 
molecules will behave as if they feel the full collective 
coupling rtfi, while for others, they will only show the 
change induced by the single-molecule coupling g. 

We thus check another observable, and ask whether the 
ground state will show correlated nuclear motion between 
distant molecules, as observed in the dark state surface. 
This can again be quantified using the expansion of the 
PES in Eq. (11). Doing so reveals a very small coupling 
parameter f3 that to lowest order is proportional to 5^/Wc 
(close to zero detuning, ojc ^ Ee{R) — Eg{R)). This corre¬ 
sponds to an even higher-order correction than the already 
small energy or bond-length shifts. Furthermore, like the 
bond-length shift, it depends on the single-molecule cou¬ 
pling instead of the collective coupling strength. We can 
thus conclude that in contrast to the excited states, the 
ground-state nuclear motion of the molecules does not 
become correlated even in the ultrastrong coupling limit. 


splitting becomes a significant fraction of the transition 
energy. Using our numerical calculations and a simple 
analytical model, we showed that while the ground-state 
energy shift is affected by the collective Rabi frequency 
(which is enhanced by VN for N molecules), other prop¬ 
erties such as the ground-state bond length depend on the 
single-molecule coupling strength and are not significantly 
affected for experimentally realistic parameters. 

Our results also lay the groundwork for a further in- 
depth exploration of the modification of molecular prop¬ 
erties under strong coupling. In particular, they provide a 
simple picture that can be used to understand the modifi¬ 
cation of chemical reactions, e.g., by lowering a potential 
barrier along a reaction coordinate. There are also a 
number of obvious extensions of the simple model pre¬ 
sented here that will be explored in the future. These 
include more realistic models of organic molecules using 
more degrees of freedom (for example, employing the PES 
obtained using quantum chemistry packages), and the 
inclusion of incoherent processes such as decay and deco¬ 
herence. We note that there has been some recent progress 
on combining QED with density functional theory [46, 47], 
which could provide complementary information to the 
model presented here. 
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V. CONCLUSIONS & OUTLOOK 

We have presented a simple model that offers a new 
perspective on strong coupling with organic molecules. 
We have shown under which conditions the molecular 
properties under strong coupling can be understood by 
the modification of the potential energy surfaces deter¬ 
mining nuclear dynamics under the Born-Oppenheimer 
approximation. In particular, we found that in many cases 
of experimental interest where the Rabi splitting is large, 
the BOA is applicable and provides an intuitive picture 
of the strongly coupled dynamics. However, we have also 
shown that for molecules with fast vibrational modes and 
large phonon-exciton couplings, the BOA can break down 
and a more involved picture is necessary. We furthermore 
demonstrated that the non-BO coupling terms between 
PES in this case are dominantly due to the change of char¬ 
acter between light and matter excitations which can be 
obtained from simple few-level models without requiring 
knowledge of the electronic wavefunctions. 

In addition, we have shown that under collective strong 
coupling involving more than one molecule, the nuclear dy¬ 
namics of the molecules in electronic “dark states” that are 
only weakly coupled to the photonic mode are nonetheless 
affected by the formation of strong coupling. In particu¬ 
lar, we find that the dark state PES describes coupling 
between the nuclear degrees of freedom of the different 
molecules. 

Finally, we investigated the change of the ground state 
properties under ultrastrong coupling, where the Rabi 
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Appendix A: Model for non-Born-Oppenheimer 
corrections 


In this appendix, we derive an analytical model for the 
non-Born-Oppenheimer corrections under molecular 
strong coupling, for a single molecule. We treat the two 
PES in the single-excitation subspace, Eg{R) + ujc and 
Ee{R), coupled by the term gHegiR). This leads to a 2 x 2 
BO Hamiltonian of the form 


HiR) = 


f Eg{R) + iOc gfiegiR)\ 

V 9tJ-eg{R) Ee{R) ) ’ 


(Al) 


which can be easily diagonalized to obtain polariton eigen¬ 
states |-|-) = cos^l^l) -|-sin0|eO) and |—) = sin^jgil) — 
cos0|eO), where |an) is short for \(j)a[x]R),n), and 


tan 29 = 


2h{R) 
6E{R) ’ 


(A2) 


where we defined 6E{R) = Eg{R) + cUc — Ee{R) and 
h{R) = ggieg{R). Using Aavg(i?) = W ^ 
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eigenenergies are given by 

E±{R) = ± ^ V4/i2(i?) + SEiR)^ . (A3) 


We can now evaluate the non-Born-Oppenheimer cou¬ 
pling terms = (fc|T„|A:') -|- (fc|^|fc')P within this 

model, using a series of approximations to obtain simple 
analytical results. First, we linearize 6E{R) « ao{R — Rc) 
around the point of intersection between the two PES, 
where Eg{Rc) + ijJc = Ef.{Rc)- Second, in the spirit of the 
Franck-Condon approximation, we assume that the dipole 
coupling is constant over the range of relevant P-values, 
and set h{R) = /iq. Following the same idea, we addi¬ 
tionally assume that the electronic wave functions do not 
change significantly with R, i.e., ^\(l)a{x;R)) « 0. This 
implies that the change in the polaritonic eigenfunctions 
|±) close to the avoided crossing at Rc is mostly due to 
the switchover between the two uncoupled surfaces, i.e., 
the change in 9{R), not because of an intrinsic change of 
electronic state with R. With these approximations, the 
correction terms become 


{-\P\+) 

{-\P^\+) 

(±|P2|±) 


—mo ho 

4hl + al{R-Rc)^’ 

2aoho(P - Rc) 

{4hl + 4{R-Rc)^)^’ 

__ 

(4h2 + a2(P-P,)2)2> 


(A4a) 

(A4b) 

(A4c) 


with the diagonal terms (±|P|±) identically zero. Note 
that diagonal terms only correspond to energy shifts and 
do not induce additional transitions [30]. The non-Born- 
Oppenheimer coupling between the polariton surfaces 
has a Lorentzian shape around the avoided crossing, and 
as expected only becomes non-negligible close to it. As 
shown in Fig. 8, the non-Born-Oppenheimer corrections 
obtained from this simple model agree almost perfectly 
with those obtained from the full numerical calculation for 
our anthracene-like model molecule. The only molecule- 
specific information entering the model are the PES of 
the uncoupled molecule and the dipole moment between 
the coupled surfaces. Specifically, the electronic wave 
functions are never used here, and their derivative as a 
function of the nuclear coordinates is not required. This 
implies that this model could be used to obtain accu¬ 
rate non-BO corrections that describe the transitions 
between potential surfaces even when the full electronic 
wave functions of a molecule are not available (e.g., in 
DFT calculations). The dynamics of the molecule could 
thus be fully recovered within a potential energy surface 
picture even when the BOA per se is not applicable. 


We now exploit this model to derive a condition for 
when the BOA becomes applicable, i.e., when the non-BO 
terms become negligible. We approximate the bare molec¬ 
ular potential energy surfaces as two harmonic oscillators 
with the same vibrational frequency ujvib, but with an 

400 

200 

3 

A 0 
O 

-200 

-400 



FIG. 8. Non-Born-Oppenheimer correction terms coupling 
the “lower polariton” and “upper polariton” PES for a sin¬ 
gle anthracene-like model molecule for a coupling strength of 
g = 0.002 a.u.. Solid colored lines: results from a full numer¬ 
ical calculation. Dashed black lines: results from the model 
Eq. (A4). Note that while all results are given in atomic units, 
the units of the P and terms are not identical, and thus 
not directly comparable. 


offset in energy and equilibrium position. 


Eg{R) « 

. r.2 

^ 2 

(A5) 

EciR) s 

i (i? Ai?)2 -b AE, 

(A6) 


where without loss of generality, we have chosen the origin 
in R and E at the minimum of Eg{R). Note that this 
model exactly results from the common approximation of 
linear coupling between a single electronic excitation and 
a bosonic vibrational mode [48, 49[. Within this model, 
SE{R) = Eg{R)+u}c — Ee{R) = ao{R — Rc) is already ex¬ 
actly linear, i.e., the linearization of the energy difference 
performed above is not an approximation. The constants 
are given by oq = and Rc = ^ + The 

maximum value of K+lj^l”)!) reached at i? = Rc, is 
given by ARujPf^/(Sho). Comparing this with the energy 
splitting at that point, E+{Rc) — E-{Rc) = 2ho, gives the 
condition that ARui^i,/(IGh^) must be small compared 
to the momentum of the respective nuclear wavefunc- 
tion (due to the additional P operating on the nuclear 
wave function). The off-diagonal terms from (—|.P^|-|-) 
reach a maximum value (again relative to the detuning) 
of MAR‘^uj^J{25V5hl) at i? = i?c + ho/{MARuU. 
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